Time - Stepping and Preserving Orthonormality

1 I n t r o d u c t i o n . This work appl ies some l inear a lgebra ideas in an o rd ina ry differential equa t ion (ODE) context . We begin by summar i s ing the a p p r o p r i a t e l inear a lgebra , and then we in t roduce the O D E problem. *Received July 1995. Revised July 1996. tThis work was supported by Engineering and Physical Sciences Research Council grants GR/H94634 and GR/KS0228. TIME-STEPPING AND PRESERVING ORTHONORMALITY 25 Suppose that A E ~:~m• with m > p, has full rank. Then A has a unique polar decomposition A = UH, U E ~t m• u T u = I, H E F~ pxp symmetric positive definite, and a unique QR decomposition A = QR, Q E F~ mxp, QTQ = I, R E ~:~pXp upper triangular with positive diagonal. We refer to U as the orthonormal polar factor and to Q as the orthonormal QR factor. Now, consider approximating A by an orthonormal matrix, measuring perturbations in the Frobenius norm, defined by ]]AIIF = trace(ATA) 1/2. We let 7(A) denote the corresponding departure from orthonormality, that is, (1.1) 7(A) := min{ JJEJJF : E e ~:[mxp, (A + E)T(A + E) = I }. It is known (see, for example, [11]) that 7(A) = II U AIIF; in other words, the orthonormal polar factor is a best orthonormal approximation. A useful order-of-magnitude estimate of 7(A) can be found using the inequalities (1.2) JIATA -IIIF < < I IAT A IIIF [[AI[2 + 1 from [11]. It is possible to construct the orthonormal polar factor U from a singular value decomposition of A, but for the application described here a more efficient approach is to use a matrix iteration. If m = p and A is nonsingular, then the sequence generated by Y0 = A, y T I//+1 = (Y/+ i ), i = 0 , 1 , 2 . . . converges quadratically to U [10]. We refer to this as the Newton iteration, since it can be obtained by applying Newton's Method to y T y __ I. For m > p, if A has full rank then the Newton iteration can be applied after an initial QR decomposition: with A -QR, if R has the polar decomposition R = URHR then A has the polar decomposition A = (QUR)HR. An alternative iteration that works for any m _> p is the Schulz iteration [13]